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For a population with any given number of types, we construct a new multivariate Moran pro- 
cess with frequency-dependent selection and establish, analytically, a correspondence to equilibrium 
Lotka-Volterra phenomenology. This correspondence, on the one hand, allows us to infer the phe- 
nomenology of our Moran process based on much simpler Lokta-Volterra phenomenology, and on 
the other, allows us to study Lotka-Volterra dynamics within the finite populations of a Moran 
process. Applications to community ecology, population genetics, and evolutionary game theory are 
discussed. 



The Moran process, originally formulated in the con- 
text of population genetics [1^ , has been applied to a wide 
range of systems in the biological sciences and statisti- 
cal physics [2-6J. Multivariate Moran models prescribe 
mean-field stochastic dynamics for birth and death in a 
population with a fixed number of individuals and any 
given number of types [5^ J . The deterministic limit 
yields replicator equations that model the evolution of 
type frequencies in large populations [3^9] . For either the 
multivariate Moran process or the corresponding replica- 
tor equations, types might be alleles, as in the population 
genetics of a single locus [D HOI [11]; species, as in the neu- 
tral theory of community ecology [l2l[T3]; or strategies, 
as in evolutionary game theory p^HTT] . Here, we build on 
a previous study of frequency- independent selection [T8] 
to construct, for any given number of types, a new multi- 
variate Moran model with frequency- dependent selection 
and demonstrate, analytically, that the corresponding 
replicator equations exhibit equilibrium Lotka-Volterra 
phenomenology. The results presented here are distinct 
from the well-known equivalence of S replicator equations 
and S — 1 Lotka-Volterra equations [2l[T9]. Furthermore, 
our choice of frequency-dependent selection generates S 
replicator equations with the equilibrium phenomenol- 
ogy of S Lotka-Volterra equations. "Lotka-Volterra equa- 
tions" in this context does not refer to the special case 
of neutrally stable, two-species, predator-prey dynamics, 
but rather to the general case of mean-field determinis- 
tic dynamics for a population with any number of well- 
mixed types and linearized density-dependent growth 
rates [20 . These generalized Lotka-Volterra equations 
are frequently used to model competitive dynamics in 
ecological communities [21]. Indeed, the two-species case 
may be the simplest dynamical system to exhibit the 
niche mechanism for stabilizing coexistence, in which (i) 
demographic rates must vary among species, (ii) abun- 
dances for each species must increase when rare, and (iii) 
m^raspecific competition must exceed interspecific com- 
petition [22j. By embedding Lotka-Volterra dynamics 
within a multivariate Moran process, we obtain a theoret- 
ical framework in which equilibrium Moran phenomenol- 
ogy may be inferred from much simpler Lotka-Volterra 
phenomenology, and Lotka-Volterra dynamics may be 



studied within the finite populations of a Moran process. 
This framework offers new insights on community ecol- 
ogy, population genetics, and evolutionary game theory. 

In a multivariate Moran process for S types and N 
individuals, the allowed states are vectors of nonnegative 
integers, ft = (ni, . . . , n^), such that < < for 
each i and J^f^^ rii = N. The allowed transitions are a 
single death event immediately followed by a single birth 
event in order to maintain a total of N individuals. In 
the absence of selection, i.e. the neutral limit, per capita 
rates of birth and death are equivalent. If we ignore 
mutation and migration, the neutral transition rate for 
a death event in Type i immediately followed by a birth 
event in Type j is simply 
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where only A^ — 1 individuals are present after the death 
of i and prior to the birth of j. Dynamics are governed by 
a multivariate master equation that we can write as jTH] 
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where is the probability of state n, is a unit vector, 
and r is a dimensionless measure of time. The &ijft = 
e(A'-(ni + l))6(n^ -l), where e{x) = 0, for x < 0, and 
1 otherwise, eliminate transitions to non-allowed states. 
The stochastic process contains S absorbing states where 
a single species dominates such that rii = N for some i. 
The more familiar univariate Moran process is obtained 
from the marginal dynamics of Eq. |2] p!8l . 

For S types, the Lotka-Volterra equations can be writ- 
ten as 



where Xi = Xi{r) is the density of Type i and 
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is a density-dependent growth rate parameterized by r^, 
the intrinsic growth rate of Type i in the absence of all 
others, and a^j, the additive per capita impact of Type j 
on the growth rate of Type i. We focus for now on com- 
petition among types such that > and a^j > for all 
i and j. In community ecology, the QjH are referred to as 
"intraspecific competition strengths" and the a^j, with 
i ^ j, as "interspecific competition strengths". With- 
out loss of generality, we can re-scale the Lotka-Volterra 
equations such that 
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and if a coexisting fixed point, x'*, exists, with x-* > 
for all z, then 
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With these preliminaries, we first construct a new 
multivariate Moran process and then establish a cor- 
respondence to equilibrium Lotka-Volterra phenomenol- 
ogy. For the construction phase, we seek frequency- 
dependent selection coefficients similar in form to the 
density-dependent growth rates of the Lotka-Volterra 
equations but remaining nonnegative for all states and 
parameter values. Towards this end, the growth rates 
of a Ricker model [23 inspire our choice of frequency- 
dependent selection coefficients 

^in = exp - ^ «i/c ^ j . (7) 
and we expand on Eq. [l] to impose the transition rates 
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where various subtractions account for the death of Type 
i prior to the birth of Type j. In the context of popu- 
lation genetics, the Wifi are reproductive fftnesses, and 
the multivariate Moran process of Eq. [2]with the transi- 
tion rates of Eq. [S] provides a framework for multiallelic, 
frequency-dependent selection that is not limited by the 
usual assumptions of weak selection or symmetry under 
the exchange of alleles [24 . Symmetry under exchange- 
ability only obtains in the neutral limit where Eq. [S] re- 

1/S for all i and 
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duces to Eq. [T] for 
3- 

The deterministic limit of Eq. [2j with the transition 
rates of Eq. [S) is given by 



dr 
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where r^' = r/N. A simpliffcation is obtained in the case 
of weak asymmetries where 
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for every z, j, and n. "Weak asymmetries" in our multi- 
variate Moran process are similar to "weak selection" in a 
univariate Moran process [H [11] , but due to the weighted 
sums over Uk that appear in Eq. [Toj our asymmetry re- 
quirements depend on population size, TV, as highlighted 
below in the discussion of Fig. [2] If we assume suffi- 
ciently weak asymmetries, and if we transform variables 
from densities, n^, to frequencies, pi = rii/N, Eq. [9]can 
be approximated by the replicator equation 
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and Wi{p) is defined by Eq. [7|upon substituting pkN for 
Uk' The mean replicator fitness is unity at all times, and 
therefore, never decreases - a standard requirement of 
replicator dynamics [2 . 

To establish a correspondence to equilibrium Lotka- 
Volterra phenomenology, we note that a coexisting fixed 
point for the Lotka-Volterra system of Eq. [3] is also a 
coexisting fixed point for the replicator system of Eq. [TT] 
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But in addition to this identity, we find a correspondence 
in stability: if the matrix with elements given by a[j -\- 
is positive definite, then the coexisting fixed point 
is globally stable in the Lotka-Volterra system and at 
least locally stable in the replicator system. The proof 
for any given number of types employs the well-known 
Lyapunov function for Lotka-Volterra dynamics [25l [26] 
(see Supplemental Material). In the special case of = 2 
competition, if r[a2i < ^2^11 and r2a'i2 < ^[^22^ then the 
replicator and Lotka-Volterra equations share a stable 
coexisting fixed point and each type can invade when 
rare. The latter inequalities, which imply that a'i2<^2i < 
^11^22 5 quantify the three requirements of niche theory 
that were highlighted in the introduction. 

To illustrate how Lotka-Volterra phenomenology can 
facilitate inference on our multivariate Moran process, 
we consider a two- type population where Eq.[2]can be re- 
written as a univariate birth-death process [18l [27] with 
rates of gain and loss given by 
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FIG. 1: Plots of the integrated marginal distribution for Type 1 in a two- type population after conditioning against extinction 
and dominance, Pn^i^) = ^ni(T)/(l — Po{r) — PN{r)). For all plots, N = 100, r[ — r2 — 0.5, and the initial condition 
is Pn/2{^ = 0) = 1. The panels correspond to scenarios of (a) neutrality {an = = = ^12 = 0-5), (b) interspecific 
exceeding intraspecific competition {a2i,ai2 = 0.52 > aii,a22 = 0.48), and (c) intraspecific exceeding interspecific competition 
(a2i,ai2 = 0.48 < aii,a22 = 0.52). At long times, the conditional probabilities approach quasi-stationary distributions 
with shapes anticipated by the stability of the coexisting Lotka-Volterra fixed point which is (a) neutrally stable (half-filled 
semicircle), (b) unstable (empty semicircle), and (c) stable (filled semicircle). 



and n = {ni^N — ni). Starting from a known ini- 
tial abundance, Fig. [l] integrates the marginal distribu- 
tion of Type 1 conditioned against extinction and dom- 
inance. We consider three cases: a) the neutral limit, 
b) destabilized coexistence due to interspecific exceeding 
intraspecific competition, and c) stabilized coexistence 
due to intraspecific exceeding interspecific competition. 
The quasi-stationary distribution emerging at long times 
is (a) fiat, (b) peaked at extinction and dominance, or 
(c) peaked at coexistence. In the corresponding Lotka- 
Volterra equations, the stability of the coexisting fixed 
point at rii/N = 0.5 allows us to infer the shape of the 
quasi-stationary distribution with (a) neutral stability 
linked to a flat distribution, (b) instability linked to a 
local minimum, and (c) stability linked to a local maxi- 
mum [27]. As N becomes large, fluctuations become rare 
and the peaks in Figs. [TJd and [T]^ approach delta func- 
tions. If we interpret Fig. [l] in the context of diploid 
population genetics, we find that interspecific exceed- 
ing intraspecific competition yields a scenario of "un- 
derdominance" (Fig. [T]3), in which selection favors ho- 
mozygous over heterozygous populations, and intraspe- 
cific exceeding interspecific competition yields a scenario 
of "overdominance" (Fig. [l]^), in which selection favors 
heterozygosity [10]. The Supplemental Material provides 
predator-prey and multi-type examples to further illus- 
trate the use of Lotka-Volterra dynamics to anticipate 
the phenomenology of our multivariate Moran process. 
For the mult i- type example, Moran dynamics are simu- 
lated with the Gillespie algorithm [28] as implemented in 
the Python package StoMPy 

But we also need to consider cases where the assump- 
tions of Eq. 
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FIG. 2: A plot of fixation rate, R, as a function of the number 
of individuals, in a two-type population where ri = r2 = 
0.5, ai2 = 0.52, and asi = 0.52 (light grey), asi = 0.56 
(dark grey), or a2i = 0.58 (black), with a[i = 1 — a2i and 
0^22 — ^ ~ Ob\2- Type 1 is always excluded as A/" ^ 00, but for 
sufficiently small A^, the a2\ > ai2 asymmetry is sufficiently 
strong that frequency-dependent selection yields fixation rates 
above the neutral limit (dashed line). 

to study Lotka-Volterra dynamics within finite popula- 
tions. As an example, consider fixation rates, or the ex- 
pected rate at which individuals of a given type grow to 
dominate a population given a fixed rate of repeated in- 
troduction. In a two-type population, the fixation rate 
of Type 1 is given by [16] 
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break down, and the correspondence 
between Moran process and Lotka-Volterra system no 
longer holds. Fortunately, in this regime, the Moran pro- 
cess remains valid for strong asymmetries and allows us 
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in our nondimensionalized time units. For large pop- 
ulations that satisfy the weak asymmetry requirements 
of Eq. 10, fixation rates are vanishingly small when 
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"^2^12 > ^i<^22 5 a^^d as expected from Lotka-Volterra dy- 
namics, Type 1 cannot invade Type 2. However, in 
the stochastic dynamics of smah populations with strong 
asymmetries, the introduction of a single individual of 
Type 1 can disproportionately impact the fitness of Type 
2 if a'21 > a'i2' For three different values of a'21 — 
Fig. [2] plots fixation rates against population sizes in a 
scenario where r2a'i2 > ^1^22- When a2i — a'12 > 0, 
fixation rates lie above the neutral limit of R = 1 for 
sufficiently small populations. Therefore, in our multi- 
variate Moran process, strong asymmetries in finite pop- 
ulations can favor invasion even when the corresponding 
Lokta-Volterra system for an infinite population predicts 
exclusion. 

We conclude by discussing two applications of our 
framework for merging Moran and Lotka-Volterra dy- 
namics. In community ecology, the past decade has 
witnessed an enormous interest in predicting large- 
scale abundance distributions based on neutral models 
for individual interactions. The foundational work of 
Caswell [12 and Hubbell [13] applied the neutral Moran 
models of population genetics to competitive ecological 
communities, such as the canopy trees in a tropical for- 
est where light-gaps are quickly filled and the total num- 
ber of individuals may be approximated as a constant. 
For these models, types correspond to species and ge- 
netic variation within each species is ignored. Efforts to 
determine the relative importance of deterministic niche 
dynamics and stochastic neutral dynamics in structuring 
ecological communities have been limited by the lack of 
a common analytical framework. Fortunately, by adding 
mutation or migration to the transition rates of our mul- 
tivariate Moran process, where selection occurs on the 
level of species rather than alleles in this context, we can 
unify Hubbell's formulation of neutral theory [T3^ with 
niche theory [22 . In particular, with the addition of a 
small mutation rate, the deterministic limit of our multi- 
variate Moran process yields a replicator-mutator equa- 
tion [161 Hi] exhibiting equilibrium Lotka-Volterra phe- 



nomenology up to small corrections. This non- neutral 
framework, in which niche stabilization can delay ex- 
tinction, offers a simple solution to Hubbell's problem 
of short species lifetimes [T3] . 

Merging Moran and Lotka-Volterra dynamics is also 
relevant to evolutionary game theory. Indeed, Fig [2] is 
reminiscent of results in Nowak et al. [30 for a univariate 
Moran process with frequency-dependent selection that 
models evolutionary games in finite populations. Their 
definition of reproductive fitness can be mapped onto 
a first-order expansion of Eq. [7] if we now assume that 
< and aij < for all i and j. In this regime, 
the equilibrium phenomenology of our re-scaled Lotka- 
Volterra and replicator equations matches the usual ex- 
pectations for two-strategy games, such as the Prisoner's 
Dilemma fSV and Repeated Prisoner's Dilemma [301132]. 
given ri = r2 and payoff matrix 

Strategy 1 Strategy 2 

strategy 1/ -an -a,2 \ (^g^ 

Strategy 2 \^ -a2i -a22 J 

Our multivariate Moran process with frequency- 
dependent selection provides a general framework for 
modeling the stochastic dynamics of evolutionary games 
with any given number of strategies in a finite popula- 
tion. 
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Supplemental Material 

A correspondence in stability between the Lotka-Volterra and 
replicator systems 

For a system of S Lotka-Volterra equations, as given by 
Eq. 3 of the main text, the well-known Lyapunov function 
can be written as 



V{x') = Ili<- <* - log^i/^D , (SM.l) 

where x[ > Q for all i, x'* is a unique global minimum, 
and the time evolution is [H [2] 



dV{x') 
dT 



1 ^ 



i=l 



(SM.2) 

Under replicator dynamics, the right-hand-side of 
Eq. SM.2, with x[ Pi, is the leading term in an ex- 
pansion of dV{p)/dT" about . Therefore, x"" is glob- 
ally stable, and ^ is at least locally stable, if a[j + a'^^ is 
positive definite. 



A predator-prey example 

The main text restricts attention to Lotka-Volterra dy- 
namics in which all the and a^j share the same sign, but 
the rescaling of Eq. 5 can be applied more generally. Of 
course, the rescaling of Eq. 5 fails when parameters are 
tuned such that r/^ = or a^i = for some i, 

but this problem occurs with vanishing probability when 



and CLij are empirically estimated. Fig.lSM.llillustrates 



our ability, in a predator-prey system, to infer the phe- 
nomenology of our multivariate Moran process from the 
corresponding Lotka-Volterra equations. A stable fixed 
point of the Lotka-Volterra system at ni/N = 0.5 antic- 
ipates the mean value of the prey's marginal stationary 
distribution. 



A multi-type example 

The main text establishes an analytical correspondence 
between our multivariate Moran process and equilibrium 
Lotka-Volterra phenomenology for any given number of 
types. Here, we provide a multi-type example in which 
the equilibrium phenomenology of our Moran process is 
anticipated by the corresponding Lotka-Volterra dynam- 
ics. We construct a consumer-resource system with S = 5 
types and J = 1000 individuals where the fifth type (con- 
sumer) benefits from the first four types (resources) . Pa- 
rameters are chosen such that, after the re-scalings of 



Predator — Prey Dynamics 
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FIG. SM.l: The integrated marginal distribution for the 
prey in a predator-prey population after conditioning against 
extinction and dominance, P^^ (r) — Pni(T)/(l — Po{r) — 
PAr(r)), for system size N = 100, initial condition P^/2(^ = 
0) = 1, intrinsic growth rates r[ = 1.25, r2 = —0.25, and 
interaction strengths a[i = 2.00, = —1.00, = ^12 — 
0.50. At long times, the conditional probability approaches 
a quasi-stationary distribution with a local maximum antici- 
pated by the stable fixed point of the Lotka-Volterra dynam- 



Eq. 5, the interaction matrix, with elements aj -, is 



/ 1.333 
0.222 
0.222 
0.222 

V -1.000 



0.211 
1.053 
0.211 
0.211 
-0.684 



0.222 
0.222 
0.889 
0.222 
-0.556 



0.222 
0.222 
0.222 
0.667 
-0.333 



0.022 \ 
0.022 
0.022 
0.022 
0.914 J 



the vector of intrinsic growth rates, with elements r^, is 

/ = (0.300, 0.300, 0.300, 0.300, -0.200), 

and all numerical values have been rounded to the near- 
est one-thousandth. Trajectories of the re-scaled Lotka- 



Volterra system are plotted in Fig. |SM.2^ with filled 
semicircles indicating equilibrium abundances that are 



globally stable due to a' 



being positive definite. 



The dynamics of our Moran process were simulated with 
the Gillespie algorithm [3 as implemented in the Python 
package StoMPy [4 . We ran 2000 stochastic trajectories 
and calculated statistics over time bins of width r"=l. 
Fig. [SM]2]d plots means and standard deviations for in- 
tegrated marginal distributions conditioned against ex- 
tinction or dominance by any single type. Clearly, at 
long-times, mean abundances of marginal stationary dis- 
tributions in the multivariate Moran process converge 
to stable equilibrium abundances of the corresponding 
Lotka-Volterra dynamics. Even over short and interme- 
diate time scales, the dynamics are remarkably similar. 
While simulations of our multivariate Moran process be- 
come computationally intensive as the number of types 
becomes large, the Lotka-Volterra system remains rela- 
tively easy to integrate. 
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FIG. SM.2: A comparison of Lotka- Volterra dynamics with the corresponding multivariate Moran process in our multi-type 
example for S — ^ and J — 1000. All types have an initial abundance of 200 in both cases. Panel (a) displays deterministic 
trajectories in the Lotka- Volterra system. Panel (b) plots means (black lines) and standard deviations (light gray regions) 
for integrated marginal distributions in the multivariate Moran process after conditioning against extinction or dominance by 
any single type. Dark gray regions indicate overlap in the standard deviations. At long-times, mean abundances of marginal 
stationary distributions in the multivariate Moran process converge to stable equilibrium abundances of the corresponding 
Lotka- Volterra dynamics (filled semicircles). 
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